
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!
 

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Initialization File
! 
! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
!       (http://www.cs.vt.edu/~asandu/Software/KPP)
! KPP is distributed under GPL, the general public licence
!       (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
!     With important contributions from:
!        M. Damian, Villanova University, USA
!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



MODULE aqchem_Initialize

  USE aqchem_Parameters
  USE aqchem_Global
  
  IMPLICIT NONE

CONTAINS

!kf Moved functions and rate coefficients that only need to be calculated once
!kf per call to cloud chem from aqchem_Rates to aqchem_Initialize 



      REAL( kind=dp )FUNCTION KMTF ( ACCOM, DG, MW )

            IMPLICIT NONE
 
            REAL( kind=dp ) ACCOM    ! accommodation coefficient (unitless)
            REAL( kind=dp ) DG       ! gas molecular diffusion coef (m2/s)
            REAL( kind=dp ) MW       ! molecular weight (g/mol)
            REAL( kind=dp ) KMT, RHO1, RAD, R, PI, V
 
            RHO1 = 1.0D0             ! density of H2O (g/mL)
            RAD = DDIAM * 0.5D0      ! droplet radius (m)
            R = 8.3145D0             ! ideal gas constant (J/mol-K)
            PI = 3.1415926536
     
            V = SQRT( 8 * R * TEMP_KPP * 1000.0D0 / PI / MW ) ! m/s
            
            KMT = ( RAD*RAD ) / ( 3 * DG )            
            KMT = KMT + ( 4 * RAD / ( 3 * V * ACCOM ) )
            KMT = 1.0D0 / KMT  ! mass transfer coefficient from Schwartz, 1986
     
            KMTF = KMT * LWC_KPP / 1000.D0 ! implied division by RHO1 (=1)
 
            RETURN
 
      END FUNCTION KMTF

      REAL( kind=dp )FUNCTION KMTB ( HLCONST, ACCOM, DG, MW )

            IMPLICIT NONE
 
            REAL( kind=dp ) ACCOM    ! accommodation coefficient (unitless)
            REAL( kind=dp ) DG       ! gas molecular diffusion coef (m2/s)
            REAL( kind=dp ) MW       ! molecular weight (g/mol)
            REAL( kind=dp ) HLCONST, RHO1, RAD, R, PI, KMT, V, R2
 
            RAD = DDIAM * 0.5D0      ! droplet radius (m)
            R = 8.3145D0             ! ideal gas constant (J/mol-K)
            R2 = 0.08206D0           ! ideal gas constant (L-atm/mol-K)  
                                     ! (= R/101325)
            PI = 3.1415926536
    
            V = SQRT( 8 * R * TEMP_KPP * 1000.0D0 / PI / MW ) ! m/s
     
            KMT = ( RAD * RAD ) / ( 3 * DG )
            KMT = KMT + ( 4 * RAD / ( 3 * V * ACCOM ) )
            KMT = 1.0D0 / KMT  ! mass transfer coefficient from Schwartz, 1986
            
            KMTB = KMT / ( R2 * TEMP_KPP * HLCONST )
 
            RETURN
 
      END FUNCTION KMTB
 
      REAL( kind=dp )FUNCTION HYDF ( KF, DH )

            IMPLICIT NONE
       
            REAL( kind=dp ) DH, KF
     
            HYDF = KF * EXP( DH * DELINVT )
            HYDF = HYDF * PHI2
        
            RETURN
     
      END FUNCTION HYDF 
        
      REAL( kind=dp )FUNCTION ORG ( KORG )

            IMPLICIT NONE
       
            REAL( kind=dp ) KORG
         
            ORG = KORG
            ORG = ORG * PHI2
         
            RETURN
     
      END FUNCTION ORG        
      
      REAL( kind=dp )FUNCTION KPHOT ( MAXC, JVAL )

            IMPLICIT NONE
 
            REAL( kind=dp ) MAXC ! photolysis rate constant max
            REAL( kind=dp ) JVAL ! photolysis rate calculated externally
 
            IF ( JVAL .LT. 0.0) THEN
               KPHOT = MAXC
               IF (CZEN_KPP .GT. 0.001_dp) THEN
                  KPHOT = PHOTO*2.718281828_dp*(MAXC)*EXP(-1/CZEN_KPP)
               ELSE 
                  KPHOT = 0.0_dp
               END IF  
            ELSE
               KPHOT = JVAL
            END IF 
 
            RETURN
 
      END FUNCTION KPHOT                  

! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! 
! Initialize - function to initialize concentrations
!   Arguments :
! 
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

! SUBROUTINE Initialize ( )

      SUBROUTINE Initialize ( TEMP2, PRES_PA, TAUCLD, PRCRATE,  &
                       WCAVG, WTAVG, AIRM, ALFA0, ALFA3, GAS, &
                       AEROSOL, CTHK1, DARK, COSZEN,          &
                       SOIL_FE_FAC, CORS_FE_FAC, SOIL_MN_FAC, &
                       CORS_MN_FAC, SEAS_NA_FAC, SOIL_NA_FAC, &
                       CORS_NA_FAC, SEAS_MG_FAC, SOIL_MG_FAC, &
                       CORS_MG_FAC, SEAS_CA_FAC, SOIL_CA_FAC, &
                       CORS_CA_FAC, SEAS_K_FAC, SOIL_K_FAC, &
                       CORS_K_FAC )

      USE aqchem_Global
      USE AQ_DATA
      USE UTILIO_DEFN  

      INTEGER :: i
      REAL(kind=dp) :: x
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

!...........Arguments:


      REAL,      INTENT( IN )  :: AIRM                    ! Total air mass in cloudy layers (mol/m2)
      REAL,      INTENT( IN )  :: ALFA0                   ! Scav coef for aitken aerosol number
      REAL,      INTENT( IN )  :: ALFA3                   ! Scav coef for aitken aerosol mass
      REAL,      INTENT( IN )  :: PRCRATE                 ! Precip rate (mm/hr)
      REAL,      INTENT( IN )  :: PRES_PA                 ! Pressure (Pa)
      REAL,      INTENT( IN )  :: TAUCLD                  ! Timestep for cloud (s)
      REAL,      INTENT( IN )  :: TEMP2                   ! Temperature (K)
      REAL,      INTENT( IN )  :: WCAVG                   ! Liquid water content (kg/m3)
      REAL,      INTENT( IN )  :: WTAVG                   ! Total water content (kg/m3)
      LOGICAL,   INTENT( IN )  :: DARK                    ! DARK = TRUE is night,  DARK = FALSE is day     
      REAL( 8 ), INTENT( IN )  :: GAS    ( : )            ! Gas phase concentrations (mol/molV)
      REAL( 8 ), INTENT( IN )  :: AEROSOL( :, : )         ! Aerosol concentrations (mol/molV)
      REAL( 8 ), INTENT( IN )  :: CTHK1
      REAL,      INTENT( IN )  :: COSZEN                  ! Cosine solar zenith angle 
      REAL( 8 ), INTENT( IN )  :: SOIL_FE_FAC             ! Fe molar fraction of ASOIL
      REAL( 8 ), INTENT( IN )  :: CORS_FE_FAC             ! Fe molar fraction of ACORS
      REAL( 8 ), INTENT( IN )  :: SOIL_MN_FAC             ! Mn molar fraction of ASOIL
      REAL( 8 ), INTENT( IN )  :: CORS_MN_FAC             ! Fe molar fraction of ACORS
      REAL( 8 ), INTENT( IN )  :: SEAS_NA_FAC             ! Na molar fraction of ASEACAT
      REAL( 8 ), INTENT( IN )  :: SOIL_NA_FAC             ! Fe molar fraction of ASOIL
      REAL( 8 ), INTENT( IN )  :: CORS_NA_FAC             ! Fe molar fraction of ACORS
      REAL( 8 ), INTENT( IN )  :: SEAS_MG_FAC             ! Na molar fraction of ASEACAT
      REAL( 8 ), INTENT( IN )  :: SOIL_MG_FAC             ! Fe molar fraction of ASOIL
      REAL( 8 ), INTENT( IN )  :: CORS_MG_FAC             ! Fe molar fraction of ACORS
      REAL( 8 ), INTENT( IN )  :: SEAS_CA_FAC             ! Na molar fraction of ASEACAT
      REAL( 8 ), INTENT( IN )  :: SOIL_CA_FAC             ! Fe molar fraction of ASOIL
      REAL( 8 ), INTENT( IN )  :: CORS_CA_FAC             ! Fe molar fraction of ACORS
      REAL( 8 ), INTENT( IN )  :: SEAS_K_FAC              ! Na molar fraction of ASEACAT
      REAL( 8 ), INTENT( IN )  :: SOIL_K_FAC              ! Fe molar fraction of ASOIL
      REAL( 8 ), INTENT( IN )  :: CORS_K_FAC              ! Fe molar fraction of ACORS 
      
!...........Local Variables:

      REAL( 8 ) :: TWASH
      REAL( 8 ) :: FECOR 
      REAL( 8 ) :: MNCOR 
      REAL( 8 ) :: NACOR 
      REAL( 8 ) :: MGCOR 
      REAL( 8 ) :: CACOR 
      REAL( 8 ) :: KCOR 
      REAL( 8 ) :: SUMPOS
      REAL( 8 ) :: SUMNEG 
      REAL( 8 ) :: Kw 
      REAL( 8 ) :: CHGBAL 

      CHARACTER( 16 ), SAVE :: PNAME = 'AQCHEM_INIT'
      CHARACTER( 120 ) :: XMSG = ' '  ! Exit status message
 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

  CFACTOR = 1.000000e+00_dp

  x = (0.0)*CFACTOR
  DO i = 1, NVAR
    VAR(i) = x
  END DO

  x = (0.0)*CFACTOR
  DO i = 1, NFIX
    FIX(i) = x
  END DO

! constant rate coefficients
! END constant rate coefficients

! INLINED initializations


      TSTART = 0.0D0
      DT     = TAUCLD
      TEND = TSTART + DT
      RTOL = 1.0D-02
      ATOL = 1.0D-02  ! was 1.d-02
     
      RTOL( ind_L_H2O2 ) = 1.0d-3
      RTOL( ind_L_GLY ) = 1.0d-3
      RTOL( ind_L_MGLY ) = 1.0d-3
      RTOL( ind_L_IETET ) = 1.0d-3 

      DDIAM = 1.6D-5 ! Droplet diameter (meters)

      TEMP_KPP = TEMP2
      PRESS = PRES_PA / 101325.D0 ! Pressure (ATM)
      LWC_KPP = WCAVG
      LWCFRAC =  LWC_KPP * 1.D-3 ! L H2O / L AIR
       
      INV_TEMP = 1.0D0 / TEMP_KPP
      DELINVT = ( 298.D0 - TEMP_KPP ) / ( 298.D0 * TEMP_KPP )

      TWASH = WTAVG * 1000.0D0 * CTHK1 * 3600.0D0 / &
             ( 1000.0D0 * MAX( 1.0D-20, REAL( PRCRATE, 8 ) ) )
          
      IF(PRCRATE .GT. 0.0) THEN
         WETFAC_KPP = 1.D0 / TWASH
      ELSE
         WETFAC_KPP = 0.D0
      END IF
      
      ALFA3_KPP = ALFA3
      ALFA0_KPP = ALFA0

      PHI2 = 1000.0D0 / 6.022d23 / LWCFRAC
      INVPHI2 = 1 / PHI2
      
      CZEN_KPP = REAL( COSZEN, 8 )

! Fraction partitioning to FE(III) and MN(II)

      IF ( DARK ) THEN
         FE_III = 0.9D0  ! Night time, GS 01July2011
      ELSE
         FE_III = 0.1D0  ! Day time, GS 01July2011
      END IF

      MN_II = 1.0D0              

! Solubility of Fe and Mn

      FE_SOL = 0.1D0               
      MN_SOL = 0.5D0    

! Set initial dynamic concentrations based on input gas and 
! aerosol concentrations 
        
      VAR( ind_G_SO2 )   = GAS( LSO2 )
      VAR( ind_G_HNO3 )  = GAS( LHNO3 ) + 2.0D0 * GAS( LN2O5 ) ! N2O5 transf. 
                                                               ! to HNO3
      VAR( ind_G_CO2 )   = GAS( LCO2 )
      VAR( ind_G_NH3 )   = GAS( LNH3 )
      VAR( ind_G_H2O2 )  = GAS( LH2O2 )
      VAR( ind_G_O3 )    = GAS( LO3 )
      VAR( ind_G_HCOOH ) = GAS( LFOA )
      VAR( ind_G_MHP )   = GAS( LMHP )
      VAR( ind_G_PAA )   = GAS( LPAA )
      VAR( ind_G_HCL )   = GAS( LHCL )
      VAR( ind_G_GLY )   = GAS( LGLY )
      VAR( ind_G_MGLY )  = GAS( LMGLY )
!      VAR( ind_G_HO   )  = GAS( LHO )
      VAR( ind_G_IEPOX ) = GAS( LIEPOX )
      
      IF( ISPC8 .gt. 0 ) THEN
      VAR( ind_G_IMAE )  = GAS( LIMAE )
      VAR( ind_G_IHMML ) = GAS( LIHMML )
      END IF
      
      VAR( ind_G_NO2 )   = GAS( LNO2 )
      VAR( ind_G_HONO )  = GAS( LHONO )
      VAR( ind_G_HNO4 )  = GAS( LHNO4 )
      VAR( ind_G_NO3 )  = GAS( LNO3RAD )
      VAR( ind_G_CH3O2 )  = GAS( LCH3O2 )
      VAR(ind_G_GCOL) = GAS( LGCOL)
      VAR(ind_G_CCOOH) = GAS( LCCOOH )
      VAR(ind_G_HCHO) = GAS( LHCHO )
      VAR(ind_G_HO2) = GAS( LHO2 ) 
      VAR(ind_G_HCHOP) = GAS( LHCHOP )
      
      IF( MTPYRAC .gt. 0 ) THEN
      VAR(ind_G_PYRAC) = GAS( LPYRUV )  
      END IF
         
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
      VAR( ind_A_NO3AKN ) = AEROSOL( LNO3, AKN )
      VAR( ind_A_NH4AKN ) = AEROSOL( LNH4, AKN )
      VAR( ind_A_CLAKN )  = AEROSOL( LCL, AKN )
      VAR( ind_A_NAAKN )  = AEROSOL( LNA, AKN )
      VAR( ind_A_SO4AKN ) = AEROSOL( LSO4, AKN )
      VAR( ind_A_PECAKN ) = AEROSOL( LEC, AKN )
      VAR( ind_A_POAAKN ) = AEROSOL( LPOA, AKN )
      VAR( ind_A_PRIAKN ) = AEROSOL( LPRI, AKN )
     
! Instantaneous droplet activation of ACC and COR modes and dissolution of H2SO4
! Initial dynamic aqueous species represent the sum of ACC and COR mode 
! concentrations.  

      VAR( ind_L_SO4MIN2 ) = AEROSOL( LSO4, ACC ) + AEROSOL( LSO4, COR )
      VAR( ind_L_NO3MIN )  = AEROSOL( LNO3, ACC ) + AEROSOL( LNO3, COR )
      VAR( ind_L_NH4PLUS ) = AEROSOL( LNH4, ACC ) + AEROSOL( LNH4, COR )
      VAR( ind_L_CLMIN )   = AEROSOL( LCL, ACC ) + AEROSOL( LCL, COR )
      VAR( ind_L_PRIACC )  = AEROSOL( LPRI, ACC )
      VAR( ind_L_NAPLUS )  = AEROSOL( LNA, ACC )
      VAR( ind_L_CAPLUS2 ) = AEROSOL( LCAACC, ACC )
      VAR( ind_L_MGPLUS2 ) = AEROSOL( LMGACC, ACC )
      VAR( ind_L_KPLUS )   = AEROSOL( LKACC, ACC )
      VAR( ind_L_PECACC )  = AEROSOL( LEC, ACC )
!      VAR( ind_L_ORGC )    = AEROSOL( LORGC, ACC )
      VAR( ind_L_ORGC )    = 0.2*AEROSOL( LORGC, ACC )
      VAR( ind_L_OXLACMIN2 ) = 0.8*AEROSOL( LORGC, ACC ) / & 
                             ( 90.03 / 177. ) ! Assume 80% of AORGC is oxalate --
			                      ! based on average results from CMAQ
			                      ! simulations where cloud generated 
					      ! org acid species were tracked                                                        
      VAR( ind_L_POAACC )  = AEROSOL( LPOA, ACC )
      
      IF( ISPC8 .gt. 0 ) THEN
      VAR( ind_L_IETET )   = AEROSOL( LIETET, ACC )
      VAR( ind_L_IEOS )    = AEROSOL( LIEOS, ACC )
      VAR( ind_L_DIMER )   = AEROSOL( LDIMER, ACC )
      VAR( ind_L_IMGA )    = AEROSOL( LIMGA, ACC )
      VAR( ind_L_IMOS )    = AEROSOL( LIMOS, ACC )
      ELSE
      VAR( ind_L_ISO3 )   = AEROSOL( LISO3, ACC )
      END IF
         
      VAR( ind_L_SO4MIN2 ) = VAR( ind_L_SO4MIN2 ) + GAS( LH2SO4 )

! Coarse crustal species from SOILCOR, ANTHCOR, SEASCOR
 
      FECOR   = SOIL_FE_FAC * AEROSOL(LSOILC,COR) + CORS_FE_FAC &
              * AEROSOL(LANTHC,COR)
      MNCOR   = SOIL_MN_FAC * AEROSOL(LSOILC,COR) + CORS_MN_FAC &
              * AEROSOL(LANTHC,COR)
      NACOR   = SEAS_NA_FAC * AEROSOL(LSEASC,COR) + SOIL_NA_FAC &
              * AEROSOL(LSOILC,COR)  &
              + CORS_NA_FAC * AEROSOL(LANTHC,COR)
      MGCOR   = SEAS_MG_FAC * AEROSOL(LSEASC,COR) + SOIL_MG_FAC &
              * AEROSOL(LSOILC,COR)  &
              + CORS_MG_FAC * AEROSOL(LANTHC,COR)
      CACOR   = SEAS_CA_FAC * AEROSOL(LSEASC,COR) + SOIL_CA_FAC &
              * AEROSOL(LSOILC,COR)  &
              + CORS_CA_FAC * AEROSOL(LANTHC,COR)
      KCOR    = SEAS_K_FAC  * AEROSOL(LSEASC,COR) + SOIL_K_FAC  &
              * AEROSOL(LSOILC,COR)  &
              + CORS_K_FAC  * AEROSOL(LANTHC,COR)
            
! Cloudwater cations 

      VAR( ind_L_NAPLUS )  = VAR( ind_L_NAPLUS ) + NACOR
      VAR( ind_L_CAPLUS2 ) = VAR( ind_L_CAPLUS2 ) + CACOR
      VAR( ind_L_MGPLUS2 ) = VAR( ind_L_MGPLUS2 ) + MGCOR
      VAR( ind_L_KPLUS )   = VAR( ind_L_KPLUS ) + KCOR            

! Fe3+ and Mn2+ for SIV oxidation
               
      VAR( ind_L_FEPLUS3 ) = FE_SOL * FE_III &
                           * ( AEROSOL( LFEACC, ACC ) + FECOR )     
      VAR( ind_L_MNPLUS2 ) = MN_SOL * MN_II &
                           * ( AEROSOL( LMNACC, ACC ) + MNCOR )
      
! Convert concententrations from input units mol/mol air --> molec/cm3 air     

      CFACTOR = PRES_PA / ( TEMP_KPP * 8.314510 ) * PHI

      DO I = 1, NVAR
         VAR( I ) = CFACTOR * VAR( I )    ! Dynamic species
      END DO

      DO I = 1, NFIX
         FIX( I ) = CFACTOR * FIX( I )    ! Fixed species
      END DO
    
!     FIX( indf_L_H2O ) = 1.D0 * INVPHI2   ! set to 1 M -- Kw includes [H2O]
      FIX( indf_L_H2O ) = 55.5D0 * INVPHI2   ! Kw in dynamic calcs does not 
                                             ! include [H2O]
!!      FIX( indf_L_HO )  = ( GAS( LHO ) * PRESS * 30.D0 * EXP( 4.5D3 * &
!!                            DELINVT ) ) * INVPHI2
!      FIX( indf_L_HO )  = INVPHI2 * ( ( GAS( LHO ) * PRESS * HOH ) /  &
!                          ( 1.d0 + HOH * WCAVG * 0.08206D0 * TEMP_KPP &
!                          / 1000.d0 ) )   !initial HO(aq) calculated based 
!                                          !on Henry's Law equil. 
!                                          !with initial/input HOg 
!                                          !representing amount of total 
!                                          !HO (gas+aq) available 
      FIX( indf_G_HO )  = GAS( LHO )*CFACTOR
      FIX( indf_L_O2 )  = 0.21 * PRESS * 1.3D-3 * EXP( 1500.D0 * DELINVT ) * INVPHI2
      
!  Calculate initial H+ and OH- from electroneutrality and Kw
!
!  Sum of positive ions + H+ = Sum of negative ions + OH-
!
!  [H+] * [OH-] = Kw 
!  [H+] = Kw / [OH-]
!
!  SUMPOS + Kw / [OH-] = SUMNEG + [OH-] -->  
!                        [OH-]**2 - (SUMPOS - SUMNEG)*[OH-] - Kw
!  Solve for [OH-] with quadratic formula and use Kw to get initial [H+]

      Kw = 1.00D-14 * EXP( -6.955D+03 * DELINVT )  ! includes conc of H2O

      SUMPOS = 2.D0 * ( VAR( ind_L_CAPLUS2)  + VAR( ind_L_MGPLUS2 ) ) &
             + VAR( ind_L_NAPLUS ) + &
               VAR( ind_L_KPLUS ) + VAR( ind_L_NH4PLUS )
      SUMNEG = 2.D0 * ( VAR( ind_L_SO4MIN2 ) + VAR(ind_L_OXLACMIN2) ) + VAR( ind_L_NO3MIN ) &
             + VAR( ind_L_CLMIN )
    
      SUMPOS = SUMPOS * PHI2
      SUMNEG = SUMNEG * PHI2
         
      CHGBAL = SUMPOS - SUMNEG
    
      VAR( ind_L_OHMIN ) = ( CHGBAL + SQRT( CHGBAL * CHGBAL &
                         + 4.D0 * Kw ) ) * 0.5D0

      IF( VAR( ind_L_OHMIN ) .gt. 0.D0 ) THEN  ! If [OH-] is (+), calculate [H+]
          VAR( ind_L_HPLUS ) = Kw / VAR( ind_L_OHMIN )
      ELSE
         XMSG = 'ERROR IN INITIAL PH CALC -- [OH-] <= 0.'
         CALL M3EXIT ( PNAME, JDATEKPP, JTIMEKPP, XMSG, XSTAT2 ) 
      END IF 
          
      VAR( ind_L_OHMIN ) = VAR( ind_L_OHMIN ) * INVPHI2   ! convert to molec/cm3
      VAR( ind_L_HPLUS ) = VAR( ind_L_HPLUS ) * INVPHI2   ! convert to molec/cm3

    

! End INLINED initializations

  RCONST(1) = ((KMTF(0.11D0,1.28D-5,64.064D0)))
  RCONST(2) = ((KMTF(0.0868D0,1.32D-5,63.013D0)))
  RCONST(3) = ((KMTF(0.00015D0,1.55D-5,44.01D0)))
  RCONST(4) = ((KMTF(0.091D0,2.3D-5,17.031D0)))
  RCONST(5) = ((KMTF(0.1532D0,1.46D-5,34.015D0)))
  RCONST(6) = ((KMTF(0.1D0,1.48D-5,47.998D0)))
  RCONST(7) = ((KMTF(0.0229D0,1.53D-5,46.025D0)))
  RCONST(8) = ((KMTF(0.006758D0,1.31D-5,48.04D0)))
  RCONST(9) = ((KMTF(0.019D0,1.02D-5,76.05D0)))
  RCONST(10) = ((KMTF(0.1158D0,1.89D-5,36.461D0)))
  RCONST(11) = ((KMTF(0.023D0,1.15D-5,58.04D0)))
  RCONST(12) = ((KMTF(0.023D0,1.15D-5,72.06D0)))
  RCONST(13) = ((KMTF(0.01D0,1.04D-5,33.01D0)))
  RCONST(14) = ((KMTF(0.0015D0,1.92D-5,46.01D0)))
  RCONST(15) = ((KMTF(0.5D0,1.3D-5,47.01D0)))
  RCONST(16) = ((KMTF(0.1D0,1.3D-5,79.01D0)))
  RCONST(17) = ((KMTB(SO2H,0.11D0,1.28D-5,64.064D0)))
  RCONST(18) = ((KMTB(HNO3H,0.0868D0,1.32D-5,63.013D0)))
  RCONST(19) = ((KMTB(CO2H,0.00015D0,1.55D-5,44.01D0)))
  RCONST(20) = ((KMTB(NH3H,0.091D0,2.3D-5,17.031D0)))
  RCONST(21) = ((KMTB(H2O2H,0.1532D0,1.46D-5,34.015D0)))
  RCONST(22) = ((KMTB(O3H,0.1D0,1.48D-5,47.998D0)))
  RCONST(23) = ((KMTB(FOAH,0.0229D0,1.53D-5,46.025D0)))
  RCONST(24) = ((KMTB(MHPH,0.006758D0,1.31D-5,48.04D0)))
  RCONST(25) = ((KMTB(PAAH,0.019D0,1.02D-5,76.05D0)))
  RCONST(26) = ((KMTB(HCLH,0.1158D0,1.89D-5,36.461D0)))
  RCONST(27) = ((KMTB(GLYH,0.023D0,1.15D-5,58.04D0)))
  RCONST(28) = ((KMTB(MGLYH,0.023D0,1.15D-5,72.06D0)))
  RCONST(29) = ((KMTB(HO2H,0.01D0,1.04D-5,33.01D0)))
  RCONST(30) = ((KMTB(NO2H,0.0015D0,1.92D-5,46.01D0)))
  RCONST(31) = ((KMTB(HONOH,0.5D0,1.3D-5,47.01D0)))
  RCONST(32) = ((KMTB(HNO4H,0.1D0,1.3D-5,79.01D0)))
  RCONST(33) = (ALFA3_KPP)
  RCONST(34) = (ALFA3_KPP)
  RCONST(35) = (ALFA3_KPP)
  RCONST(36) = (ALFA3_KPP)
  RCONST(37) = (ALFA3_KPP)
  RCONST(38) = (ALFA3_KPP)
  RCONST(39) = (ALFA3_KPP)
  RCONST(40) = (ALFA3_KPP)
  RCONST(92) = ((1.1d0))
  RCONST(96) = (WETFAC_KPP)
  RCONST(97) = (WETFAC_KPP)
  RCONST(98) = (WETFAC_KPP)
  RCONST(99) = (WETFAC_KPP)
  RCONST(100) = (WETFAC_KPP)
  RCONST(101) = (WETFAC_KPP)
  RCONST(102) = (WETFAC_KPP)
  RCONST(103) = (WETFAC_KPP)
  RCONST(104) = (WETFAC_KPP)
  RCONST(105) = (WETFAC_KPP)
  RCONST(106) = (WETFAC_KPP)
  RCONST(107) = (WETFAC_KPP)
  RCONST(108) = (WETFAC_KPP)
  RCONST(109) = (WETFAC_KPP)
  RCONST(110) = (WETFAC_KPP)
  RCONST(111) = (WETFAC_KPP)
  RCONST(112) = (WETFAC_KPP)
  RCONST(113) = (WETFAC_KPP)
  RCONST(114) = (WETFAC_KPP)
  RCONST(115) = (WETFAC_KPP)
  RCONST(116) = (WETFAC_KPP)
  RCONST(117) = (WETFAC_KPP)
  RCONST(118) = (WETFAC_KPP)
  RCONST(119) = (WETFAC_KPP)
  RCONST(120) = (WETFAC_KPP)
  RCONST(121) = (WETFAC_KPP)
  RCONST(122) = (WETFAC_KPP)
  RCONST(123) = (WETFAC_KPP)
  RCONST(124) = (WETFAC_KPP)
  RCONST(125) = (WETFAC_KPP)
  RCONST(126) = (WETFAC_KPP)
  RCONST(127) = (WETFAC_KPP)
  RCONST(128) = (WETFAC_KPP)
  RCONST(129) = (WETFAC_KPP)
  RCONST(130) = (WETFAC_KPP)
  RCONST(131) = (WETFAC_KPP)
  RCONST(132) = (WETFAC_KPP)
  RCONST(133) = (WETFAC_KPP)
  RCONST(134) = (WETFAC_KPP)
  RCONST(135) = (WETFAC_KPP)
  RCONST(136) = (WETFAC_KPP)
  RCONST(137) = (WETFAC_KPP)
  RCONST(138) = ((KMTF(0.05D0,1.53D-5,17.007D0)))
  RCONST(139) = ((KMTB(HOH,0.05D0,1.53D-5,17.007D0)))
  RCONST(140) = (WETFAC_KPP)
  RCONST(141) = ((KMTF(0.02D0,1.0D-5,118.13D0)))
  RCONST(142) = ((KMTF(0.02D0,1.0D-5,102.0D0)))
  RCONST(143) = ((KMTF(0.02D0,1.0D-5,102.0D0)))
  RCONST(144) = ((KMTB(HIEPOX,0.02D0,1.0D-5,118.13D0)))
  RCONST(145) = ((KMTB(HMAE,0.02D0,1.0D-5,102.0D0)))
  RCONST(146) = ((KMTB(HHMML,0.02D0,1.0D-5,102.0D0)))
  RCONST(158) = (WETFAC_KPP)
  RCONST(159) = (WETFAC_KPP)
  RCONST(160) = (WETFAC_KPP)
  RCONST(161) = (WETFAC_KPP)
  RCONST(162) = (WETFAC_KPP)
  RCONST(163) = (WETFAC_KPP)
  RCONST(164) = (WETFAC_KPP)
  RCONST(165) = (WETFAC_KPP)
  RCONST(166) = (WETFAC_KPP)
  RCONST(167) = ((KMTF(0.023D0,1.15D-5,60.052D0)))
  RCONST(168) = ((KMTF(0.0322D0,1.24D-5,60.052D0)))
  RCONST(169) = ((KMTF(0.02D0,1.64D-5,30.03D0)))
  RCONST(170) = ((KMTF(0.02D0,1.64D-5,30.03D0)))
  RCONST(171) = ((MTPYRAC*KMTF(0.0322D0,1.24D-5,88.06D0)))
  RCONST(172) = ((KMTB(GCOLH,0.023D0,1.15D-5,60.052D0)))
  RCONST(173) = ((KMTB(CCOOHH,0.0322D0,1.24D-5,60.052D0)))
  RCONST(174) = ((KMTB(HCHOH,0.02D0,1.64D-5,30.03D0)))
  RCONST(175) = ((KMTB(HCHOH,0.02D0,1.64D-5,30.03D0)))
  RCONST(176) = ((MTPYRAC*KMTB(PYRACH,0.0322D0,1.24D-5,88.06D0)))    
  RCONST(189) = ((HYDF(1.8D-1,4.03D+3)))
  RCONST(190) = (5.1D-3)
  RCONST(191) = ((HYDF(1.8D-1,4.03D+3)))
  RCONST(192) = (5.1D-3)
  RCONST(193) = ((KPHOT(4.6D-6,JH2O2)))
  RCONST(212) = (WETFAC_KPP)
  RCONST(213) = (WETFAC_KPP)
  RCONST(214) = (WETFAC_KPP)
  RCONST(215) = (WETFAC_KPP)
  RCONST(216) = (WETFAC_KPP)
  RCONST(217) = (WETFAC_KPP)
  RCONST(218) = (WETFAC_KPP)
  RCONST(219) = (WETFAC_KPP)
  RCONST(220) = (WETFAC_KPP)
  RCONST(221) = (WETFAC_KPP)
  RCONST(222) = (WETFAC_KPP)
  RCONST(223) = (WETFAC_KPP)
  RCONST(224) = (WETFAC_KPP)
  RCONST(225) = (WETFAC_KPP)
  RCONST(226) = (WETFAC_KPP)
  RCONST(227) = (WETFAC_KPP)
  RCONST(228) = ((KMTF(0.05D0,1.D-5,62.0D0)))
  RCONST(229) = ((KMTF(0.05D0,1.35D-5,49.0D0)))
  RCONST(230) = ((KMTB(NO3H,0.05D0,1.D-5,62.0D0)))
  RCONST(231) = ((KMTB(CH3O2H,0.05D0,1.35D-5,49.0D0)))
  RCONST(232) = ((KPHOT(5.6D-7,-1.d0)))
  RCONST(261) = (WETFAC_KPP)
  RCONST(262) = (WETFAC_KPP)
  RCONST(263) = (WETFAC_KPP)
  RCONST(264) = (WETFAC_KPP)
  RCONST(265) = (WETFAC_KPP)
  RCONST(266) = (WETFAC_KPP)
  RCONST(267) = (WETFAC_KPP)
  RCONST(268) = (WETFAC_KPP)
  RCONST(269) = (WETFAC_KPP)        
      
END SUBROUTINE Initialize

! End of Initialize function
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



END MODULE aqchem_Initialize

